Femtoscopy in hydro-inspired models with resonances * 
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Effects of the choice of the freeze-out hypersurface and resonance decays on the HBT inter- 
ferometry in relativistic heavy-ion collisions are studied in detail within a class of models with 
single freeze-out. The Monte-Carlo method, as implemented in THERMINATOR, is used to generate 
hadronic events describing production of particles from a thermalized and expanding source. All 
well-established hadronic resonances are included in the analysis as their role is crucial at large 
freeze-out temperatures. We use the two-particle method to extract the correlation functions, which 
allows us to study the Coulomb effects. We find that the pion HBT data from RHIC are fully 
compatible with the single freeze-out scenario, pointing at the shape of the freeze-out hypersurface 
where the transverse radius is decreasing with time. Results for the single-particle spectra for this 
situation are also presented. Finally, we present predictions for the kaon femtoscopy. 
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INTRODUCTION 



Femtoscopy is one of the most important and promis- 
ing techniques used in relativistic heavy-ion collisions, 
as it reveals the spacetime characteristics of the fireball 
formed in the reaction. The study of two-particle cor- 
relations of identical particles is known as the Hanbury- 
Brown-Twiss (HBT) interferometry 0,0,0,1111, while 
together with the extension to non-identicalparticles it 
has been generically termed femtoscopy 0,0j, referring 
to studies of the system at the femtometer scale. For 
a recent review of various aspects, both theoretical and 
experimental, of the field the reader is referred to the re- 
view by Lisa, Pratt, Soltz, and Wiedemann §j. There 
are several questions concerning the heavy-ion data for 
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the HBT radii. The major puzzle is the practically con- 
stant value of the HBT radii over the huge reaction en- 
ergy range ^/sjvjv = 20—200 GeV. The other surprising, 
when confronted to expectations based on earlier model 
predictions, feature found at RHIC is the proximity of 
the value of i? ut /Rside to unity, indicating a very short 
emission time of pions from the source. These challenges 
are to be faced by theoretical descriptions. As a mat- 
ter of fact, a simultaneous reproduction of all HBT radii 
poses a serious problem to models with limited paramet- 
ric freedom, as well as to hydrodynamical simulations or 
transport codes (for a recent study see pj] , where Rs and 
Rl are found to be in agreement with data, while Ro 
is predicted to be larger that the observed one). The 
above puzzles and problems led to significant revision of 
our understanding of the RHIC physics, with such im- 
portant conclusions as the absence of the latent heat in 
hadronization 3, which would lead to a much longer- 
lived fireball, or the introduction of the concept of the 
length of homogeneity 53|, which effectively reduces the 
observed radii to smaller values than the geometric size 
of the whole source. 

In the present paper we study the femtoscopy at RHIC 
in a class of hydro-inspired models, all with a single- 
freeze-out £jj . Our analysis uses THERMINATOR 12] - the 
THERMal heavy IoN generATOR for the single-freeze- 
out approach, which is a very flexible tool for studies of 
this type. The single freeze-out concept, which identi- 
fies the thermal and kinetic freeze-outs, may be viewed 
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as an approximation to a more detailed evolution, taking 
into account different time scales for various processes. 
Nevertheless, the single freeze-out complies to the explo- 
sive scenario at RHIC an d is definitely worth a de- 
tailed study in the context of femtoscopy. Moreover, the 
approach reproduces very efficiently the particle abun- 
dances, the transverse-momentum spectra, including par- 
ticles with strangeness , produces very reasonable re- 
sults for the resonance production th e balance func- 
tions in rapidity 0], the elliptic flow [17j . and the trans- 
verse energy ■ Approximate predictions of the model 
for the HBT radii were already presented in Ref. [rRll9| . 

In particular, we deal in detail with two important is- 
sues. The first one is the influence of the chosen freeze- 
out hypersurface X on the model prediction for the pio- 
nic HBT radii. We study several parametrization of X: 
the one from the original single-freeze-out model [ijj, as 
well as three hypersurfaces from the Blast- Wave model 
[2H l2l| and its extension which account for the changes 
of the shape in the t — p (time - transverse radius) space 
with the help of a single parameter a. We find signif- 
icant dependence here, especially for the i?i on g radius. 
The analysis favors the shape of the freeze-out hypersur- 
face where the transverse radius is decreasing with time. 
We also present the results for the single-particle spec- 
tra for the best freeze-out model out of the four studied, 
which turns out to be the modified Blast- Wave model 
with parameter a = —0.5. Clearly, the HBT correlation 
data place constraints on the freeze-out geometry and 
our study helps to put quantitative bounds on model pa- 
rameters. 

The second issue studied is the detailed role of hadronic 
resonances on the two-particle pion correlation. Al- 
though the basic features have been known since the 
very early days of the field 0, there have been no 
HBT studies at large temperatures (T ~ 170 MeV) tak- 
ing into account sufficiently many resonances. The inclu- 
sion of practically all resonances at such temperatures is 
important, as revealed by the studies of particle ratios 
Hi 13 EE EE EllH and the transverse-momentum 
spectra [ll|, where the resonances significantly decrease 
the inverse-slope parameter |29j. The presence of res- 
onances plays an essential role also in femtoscopy sup- 
plying the separation distributions with long exponential 
tails, thus providing non-gaussian features in the HBT 
correlation functions. For a very recent study see Ref. 

We should stress, that the models termed here 
as "Blast- Wave" have the blast-wave geometry, however 
they do include all resonances, contrary to many appli- 
cations, where no resonance feeding is present. 

Finally, we present predictions for the kaon femtoscopy. 
We find that when the pion and kaon results are plotted 
together, they comply to the ra^-scaling conjuncture. 

The paper is constructed as follows: In Sect. II we de- 
fine a class of the hydro-inspired models with resonances 
which are used in our analysis of the correlation func- 
tions. Section III is a short introduction to the HBT 
formalism. In Sect. IV we present our main results on 




FIG. 1: Various parametrization of the freeze-out hypersur- 
face. The curves show the dependence of time t on the radial 
distance p = yrl+r| at r 2 = for the four models consid- 
ered. 

the HBT radii. The effects of the resonances on the cor- 
relation function are discussed in Sect. V. Section VII 
contains our predictions concerning the kaon correlation 
functions. 



II. HYDRO-INSPIRED MODELS OF 
FREEZE-OUT 

The hydro-inspired models have become a very popu- 
lar tool to analyze the data collected in relativistic heavy- 
ion collisions The most popular model belonging 
to this class is the Blast- Wave model of Schnedermann, 
Sollfrank and Heinz [if]]. In the original formulation, 
it was designed to describe boost-invariant and cylindri- 
cally symmetric systems, hence it is best suited for de- 
scription of the midrapidity region of the central Au + 
Au collisions studied at the top RHIC energies. 

Other models of this type include the Buda-Lund 
model [3^ 01 an d the Cracow single-freeze-out model 
[ill ll 4. . A distinctive feature of the Cracow model is 
that it includes the complete set of hadronic resonances. 
Precisely this feature allowed for the uniform description 
of the chemical and thermal freeze-outs. 

In this work we consider the boost-invariant and cylin- 
drically symmetric systems and use THERMINATOR [12] to 
include the resonance effects. In this way, the Blast- 
Wave and Cracow model are treated on the same footing; 
the contributions from the resonance decays that are very 
often neglected in the studies based on the Blast-Wave 
model are now taken completely into account. The only 
(important) difference between the considered models re- 
sides in the definition of the freeze-out hypersurface; for 
the Blast- Wave model the freeze-out hypersurface is typ- 
ically defined by the condition of the constant laboratory 
time, whereas for the Cracow model the freeze-out hyper- 
surface is defined by the condition of the constant proper 
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final pion 




FIG. 2: The cascade of resonance decays: the initial resonance 
formed on the freeze-out hypersurface at the point xn with 
momentum pjv decays after proper time rjv at the point xn-i- 
We track one decay product, which decays in sequence until 
the final pion is formed at the point x 1 . 

time. In the present paper we take into consideration 
these two cases and other possible options, all illustrated 
in Fig. H The lines show the freeze-out hypersurfaces 
at r z = 0. Due to the measure 2irpdp the parts of the 
hypersurface with larger values of p are more relevant in 
the evaluation of observables. 

The boost-invariant, cylindrically symmetric freeze- 
out models are distinguished by different freeze-out 



curves in the t — p space (t is the laboratory time mea- 
sured at r z = 0, while p = yjr% + is the distance 

from the collision axis). The most popular Blast-Wave 
parametrization uses the freeze-out condition t = const. 
In our calculations we also consider parametrizations of 
the form t = t + ap, where r and a are constants. 
For a = we reproduce the standard case. By study- 
ing the cases with different values of a we may analyze 
the effect of the shape of the freeze-out hypersurface on 
different physical observables, including the HBT radii. 
The positive (negative) values of the parameter a corre- 
spond to the freeze-out curves which go up (down) in the 
Minkowski t — p space, cf. Fig.^ It is important to real- 
ize that the freeze-out curves with negative a's resemble 
the freeze-out curves obtained typically in hydrodynamic 
calculations. In the latter case the matter placed close 
to the surface of the system decouples earlier, while the 
matter inside the system decouples later, only when the 
temperature in that region drops down due the expansion 
of the system. Consequently, by studying the effect of the 
choice of the freeze-out curve on the extracted HBT radii, 
we may find which freeze-out models are favored by the 
data. 

A. Emission function 

Our starting point is the formula for the pion emis- 
sion function which takes into account sequential decays 
of the resonances 0,0, A contribution from one 
particular decay chain c, illustrated in Fig. [3 is given by 
the following equation 



S c (x uPl )=E pl j^— = J ^*B(p2,pi) J dr 2 T 2 e- T ^ J d 4 x 2 S^ (x2 + 
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J B (pn,Pn-i) J dT N T N e TnTn J dS M (x N ) p% £ (4) (^x N + P ^J^ ~ %N-lj In [pn ■ u (x N )} . (1) 



Here dS M is a 3-dimensional element of the freeze-out hy- 
persurface S, the position of a resonance that decouples 
on the freeze-out hypersurface S is denoted by xn, its 
four-momentum by p N , and its mass by tun- The func- 
tion /at is the thermal distribution function depending on 
the product of the resonance four-momentum P n and the 
local four-velocity u(xn)- The resonance decays at the 
spacetime point xn-i and produces the tracked daugh- 
ter particle with the four-momentum p N _ 1 . In the first 
step of the simulation of the decay, the proper time tn 
is generated randomly according to the exponential de- 
cay law (1/rjy) exp(— TnTn), where Tn is the resonance 
width. Then, in the second step, the position xn-i is ob- 



tained from the formula x N _ 1 = x 1 ^ + (p n /itin)tn- The 
momentum of the daughter particle p N _ 1 is determined 
purely by the available phase space, as described in Refs. 
[la Hfll | . In Eq. the phase space effects are taken into 
account in terms of the splitting functions B (pn,Pn-i) 
defined in [19]. If the daughter particle is a resonance, 
the above scheme is repeated until the final stable par- 
ticle (pion, kaon, nucleon or antinucleon) is produced. 
The spacetime position and four-momentum of the final 
particle is denoted as x\ and P \ . The cascade character 
of the process is reflected by the iterative structure of 
Eq. GJ. The complete emission function is obtained as 
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the sum over all possible decay channels c, 

S(x,p) = ^2S c (x,p). (2) 

c 

We note that all considered particles are on the mass 
shell with the energies E p — y/ m? + p 2 . 

B. Distributions of the primordial particles 



distance from the collision axis. The function arj_(C) 
conveniently parametrizes the transverse collective flow, 
w_l(C) = tanhaa(C). We recall that the longitudinal flow 
has the form v z — r z /t The variable (3 appearing 

in Eq. <PJ is the inverse temperature, f3 — 1/T, whereas 
the plus-minus sign is related to the statistics (in Eq. © 
and in the expressions below the upper sign corresponds 
to fermions, while the lower sign to bosons). With the 
help of the modified Bessel functions, Eq. <EJ may be 
rewritten in the more compact form 



The Monte-Carlo simulation of the decay process starts 
with the generation of hadronic distributions at the mo- 
ment of freeze-out. We shall refer to such distributions 
as the primordial distributions. All hadronic states ap- 
pearing in the Particle Data Tables [36] are included in 
this procedure. Since we consider a boost-invariant and 
cylindrically symmetric system, the primordial distribu- 
tions in rapidity, y, and transverse momentum, p±, are 
of the form 



dN 1 



dyd 2 p A 



1 f 



dp 



x \ m ±-r:Ki [nf3m±cosha±] Iq [nPp±smha±] 

df 1 
— P± — Kq [n/3m^coshajJ I\ [n/3p_i_sinharij > . (4) 
dQ J 



dN 



dyd 2 pA 



dZ^XN^f (p N ■ u(x N )) 



2n oo 



(2vr 



TO±cosh(a|| -y)^ -p± cos(0 - tp)^ 



x < exp [/3mj_cosh(aj_(C))cosh(a|| — y) 



-f3p±smh(a±(()) cos(0 — ip) — (3p\ ± 1 



(3) 



Here = ^Jm 2 + p 2 ± is the transverse mass, and p 
is the momentum azimuthal angle. The quantities <f>, 
arii, and C are used to parametrize the freeze-out hyper- 
surface: <f) is the azimuthal angle, tan0 = r y /r x , a\\ is 
the spacetime rapidity, a\\ = l/21n[(t + r z )/(t — r z )], 
while C parametrizes the freeze-out curve obtained as the 
projection of the freeze-out hypersurface on the r z = 
plane. The freeze-out curve is defined by the mapping 
C - * (^(C))i°(C)) relating the freeze-out time and posi- 
tion. At r z = the variable f = ^Jt 2 — r 2 coincides 

with the laboratory time, whereas p — \/r 2 +r 2 is the 



C. Cracow single- freeze-out model 

Different boost-invariant and cylindrically symmet- 
ric models may differ in the choice of the freeze-out 
curve (t(C) i p(C))- I n the Cracow single-freeze-out model 
[ill 13 Ha] the freeze-out hypersurface is specified by the 
following equations: 



f = t cosha;_L(C), p — t sinho!j_(C), t = const, 
which are equivalent to the condition 



f 2 - p 2 = t 2 



2 2 2 

r — r = t . 

x y 



(5) 



(6) 



The parameterization iJ^J implies that the freeze-out of 
the fluid elements placed farther away from the center 
happens at later times, see Fig. ^ The velocity profile 
in the Cracow model has the Hubble-like structure 



r 

v = — . 

t 



(7) 



The use of Eqs. © and Q in Eq. JHJ and the change 
of the integration variable £ (first to a± and later to p) 
leads to the distribution of the primordial particles in the 
6-dimensional space of spacetime positions and momenta 



dN 



1 



dydipp±dp±da\\d(j)pdp (2n) 3 

x < exp 



m>x \/t 2 + p 2 cosh(a II — y) — pj_pcos(0 — ip) 



p p 

0mj_ \/ H — ^cosh(a|| — y) — f3p± — cos(</> — ip) — (3p 



± 1 



(8) 



r 



D. Generalized Blast-Wave Model 

As the second option considered in our calculations we 
choose the generalized Blast- Wave parameterization 

f = t + ap, tanhax(C) = v± = const. (9) 



The parameters: t, a and v± are constants (for simplicity 
we assume that the transverse flow profile is constant). 
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For a — we obtain the standard Blast-Wave parameter- 
ization corresponding to the assumption that the freeze- 
out process happens at constant laboratory time t = r 
(in the central region where an = 0). For a > (a < ) 



the straight line defining the freeze-out in the Minkowski 
space at r z = goes upwards (downwards). Similarly to 
the previous case, the use of the parameterization ^ in 
Eq. <PJ gives the 6-dimensional density 



dN 



dydtpp±dp±da\\ dtftpdp 




(r + ap) [mj_cosh(a|| — y) — ap± cos(<p — ip)] 
(3m 



cosh(a|| — y) -j=k ± cos(0 — ip) — f3p 



± l 



(10) 



III. THE HBT FORMALISM 

A. Two-particle correlation function 

Consider the two-particle distribution expressed by the 
two-particle emission function, 



W 2 (Pi,P 2 ) 



F F dN 
Pl P2 d^p 2 



S{x 1 ,x 2 ,p 1 ,p 2 )d A x 1 d A x 2 . (11) 



The correlation function is then defined as 

W 2 (p 1 ,p 2 ) 



where 



C(Pi,P 2 ) 



W 1 (p)=E 1 



W 1 {p 1 )W 1 {p 2 ) 



dN 
d 3 p 



d 4 x S(x,p) 



(12) 



(13) 



with S(x,p) given by Eqs. and (2J. 

One may assume that the two-particle production 
probability is influenced only by the two-particle inter- 
action. In this case one neglects the many-body inter- 
actions between the produced particles as well as the 
event-wide correlations (e.g., the effects induced by the 
momentum conservation). Then, the two-particle emis- 
sion function may be expressed as the product of the 
single-particle emission functions and the squared wave- 
function of the pair. After taking into account the 
smoothness approximation we write 



C(q,k) 



J d i x 1 S(x uPl )d' l x 2 S(x 2 ,p 2 )\^(k*,r 
J d 4 XiS(xi,pi) J d 4 x 2 S{x 2 ,p 2 ) 



*M2 



(14) 

We define the momentum difference of the particles as 

q = 0o,g) = (E P1 -E P2 ,p 1 -p 2 ), (15) 

the sum of their momenta as 

P=(P ,P) = (E pl +E P2 , Pl +p 2 ), (16) 



and the average momentum of the pair as 

k = {k Q ,k) = -{E Pl +E P2lPl +p 2 ). (17) 

The generalized momentum difference is defined by the 
formula 



q = q 



P(q ■ P) 
p2 ■ 



(18) 



which in the pair rest frame (PRF) is reduced to the form 



9=(0,2fc*) 



(19) 



For the particles with equal masses we use the notation 



= 2k* 



(20) 



The space and time separations of the members of the 
pair are: r = r\ — r 2 and At = t% — t 2 . If calculated in 
PRF they are denoted as r* and At*. Both k* and r* 
appear as the arguments of the wave function in Eq. l|T4*jl 
since PRF is the most convenient reference frame for the 
representation of the wave function structure. 

In general, the HBT analysis may be performed in any 
reference frame. One determines the correlation function 
as a function of the relative-momentum components in 
the selected frame. Then the inverse widths of the corre- 
lation functions yield the size parameters of the system 
in this frame. In the present paper we use the Bertsch- 
Pratt decomposition (slj. l^ol l4(l | of the average and rela- 
tive three-momenta into three components. The long axis 
coincides with the beam axis, the out axis is determined 
by the direction of the average transverse momentum of 
the pair, denoted later by fey, and the side direction is 
perpendicular to the other two axes. Following the RHIC 
experiments we choose to perform the analysis in the lon- 
gitudinal co-moving system (LCMS), which is defined as 
a system where fci ong = 0. The HBT radii presented 
below are always obtained in the LCMS system. Later 
in this work the notation is used in which the values in 
PRF are denoted by an asterisk, while the values without 
asterisk are defined in LCMS. 
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By the definition of the Monte-Carlo method, the nu- 
merical equivalent of the integrals itTHl is the summa- 
tion over particles or pairs of particles generated by the 
Monte-Carlo procedure. The numerical calculation of the 
correlation functions is done in bins, which may be ex- 
pressed with the help of the function 

f 1 if |x| < #, 
<Ja(x) = I " (21) 

I otherwise. 

Then the correlation function may be expressed simply 
as 

C(q,k) = 

£ £ Ms -p* - ^ + P,))l*( fc >*)l 2 

£Z>A(9-P»+Pj)<*A(fc- \(Pi+Pj)) 

i 3 

(22) 

In our numerical calculations we use A = 5 MeV. 




m T - m [GeV] 



FIG. 3: Transverse-mass spectra at mid-rapidity of pions 
(open circles) and kaons (open squares) for the Blast- Wave 
model with a = -0.5, T = 165.6 MeV, \i B = 28.5 MeV, 
r = 8.55 fm, p max = 8.92 fm, and vt = 0.311 c. The data 
points (stars for pions, triangles for kaons) come from the 
STAR collaboration |41j . 



B. Wave function of the pion pair 

Various analyses of the HBT correlations use different 
approximations for the full pion wave function 4'. In 
the non-interacting system or in the interacting but non- 
relativistic case the motion of the center of mass can be 
separated and one deals with the relative motion only. 
The simplest relative wave function ignores all dynamical 
interactions and has the form 



+ e 



— ifc* i 



(23) 



where symmetrization over the two identical particles has 
been performed. Therefore, 



|*0| 2 = l + cos(2fc*r*). 



(24) 



Correlation functions calculated according to (l22l) and 
Il24t represent the ideal Bose-Einstein correlation func- 
tions. They are also very useful in the model studies, 
because they can be calculated analytically for simple 
gaussian emission functions. 

In more realistic calculations, the Coulomb interaction 
of the charged pion pairs should be taken into account, 
which may be achieved by the use of the wave function 



V2 



r F(-iTI,l,it+) 



+ e*'*' F(-i V ,l,iC) 



(25) 



where S c is the Coulomb phase shift, A c is the Coulomb 
penetration factor (sometimes called the Gamow factor) , 
£± = k*r*±k*r* = k*r* (1 ± cos 9*), T] = (k*^- 1 with a 



being the Bohr radius of the pair, and F is the confluent 
hypergeometric function. 9* is the angle between k* and 
r* . The correlation function obtained in this way can be 
compared directly to the correlation function obtained 
from the experiment. The complete wave function of 
the pion pair contains a contribution from the strong 
interaction as well. However these are small in the isospin 
1 = 2 channel and are neglected. 



C. Numerical calculation of correlation functions 

The correlation functions analyzed in this work are ob- 
tained through a numerical implementation of Eqs. l(22|l 
and ll23"j) or Eqs. l|22j) and Particles generated by 

THERMINATOR are grouped into events, as in experiment. 
In each event every charged pion is combined with every 
other pion of the same charge. For each pion pair, l^l 2 is 
calculated and added to the numerator of Eq. ll22ll in a 
bin corresponding to its q mv for 1-dimensional functions 
or to its g ut) 9sidc and qWg for the full 3-dimensional 
case. At the same time, 1 is added to the denominator 
of Eq. (52} in the corresponding bin. The resulting ratio 
yields the correlation function. 

By making a proper selection of single pions and pairs 
of pions one may study the correlation functions as func- 
tions of various variables. For instance, taking into ac- 
count the pairs of particles within a certain total momen- 
tum range only, one immediately obtains the dependence 
on hr. It is important to note that all single-particle ap- 
proaches in the studies of the correlation functions use 
transverse-momenta of single particles only, while the ex- 
perimental data are represented as functions of fey. This 
is one of the advantages of the two-particle method over 
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k T [GeV/c] k T [GeV/c] 



FIG. 4: The results for Cracow model: A and the HBT 
radii Rout, -Rsidc, and i?i ong shown as functions of the trans- 
verse momentum of the pion pair. The squares show the 
full calculation with resonances based on the method of 
Sect. IIIII down-triangles is the same without resonances, 
the up-triangles show the calculation with resonances and the 
Coulomb corrections made according to the Bowler-Sinyukov 
method, while the circles show the data of the STAR collabo- 
ration for y/SNN = 200 GeV 01 . The lines are drawn to guide 
the eye. We note that the inclusion of resonances increases the 
radii by about 1 fm. The model parameters are: T = 165.6 
MeV, fi B = 28.5 MeV, r = 10.55 fm, and p max = 7.53 fm. 



FIG. 5: Same as Fig. [21 for the Blast-Wave model with reso- 
nances, a = 0.5. The model parameters are: T = 165.6 MeV, 
fi B = 28.5 MeV, r = 9.91 fm, p max = 7.43 fm, and v T = 
0.407 c. 



D. Extraction of HBT radii 

In most of the realistic cases the integral 11411 cannot 
be performed analytically. One of the cases where the 
analytic calculation may be done corresponds to the sit- 
uation where the pion wave function is given by Eq. \'2'M 
and the single-particle emission function is a static 3- 
dimensional ellipsoid with a gaussian density profile 



the single-particle method. Another one is the possibil- 
ity of including final-state interactions, such as Coulomb 
effects. 



S{x,p) = Nexp - 



''side 



2Hi 



2m 



x long 



(26) 



Please note that this source function is static - it does not 
depend on particle momentum. In this case the integral 
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FIG. 6: Same as Fig. [3] for the Blast-Wave model with reso- 
nances, a = 0. The model parameters are: T = 165.6 MeV, 
fi B = 28.5 MeV, r = 8.17 fm, p max = 8.21 fm, and v T = 
0.341 c. 

lllH with the free wave function Il2-'ill leads to the well 
known formula 

C (k±,q out ,q s idc, giong) = 1 + Acxp [-Rl ut (k±)ql ut 

~ ^sideC^-OOsidc — Aong(^i) < ziong] ■ (27) 

The quantities i? ut, -Rside and -Ri on g, known as the "HBT 
radii" are the widths of the gaussian approximation to 
the single-particle freeze-out distribution. It is important 
to emphasize that formula l(27l is commonly used to fit 
the experimental data and to represent the results of the 
model calculations although the experimental or model 
emission functions are frequently far from gaussians. In 
context of our model this issue will be discussed in detail 
in Sect. V. 



FIG. 7: Same as Fig. El f° r the Blast- Wave model with 
resonances, a — —0.5. The model parameters are: T = 
165.6 MeV, \i B = 28.5 MeV, r = 8.55 fm, p max = 8.92 fm, 
and vt = 0.311 c. This is the model that produces best agree- 
ment out of the four models considered. With the same values 
of the parameters the model reproduces the transverse-mass 
spectra, see Fig. El 



E. Bowler- Sinyukov formalism 

The pion correlation functions for sizes typically ob- 
served in heavy-ion collisions are significantly influenced 
by the Coulomb interaction, hence the model calcula- 
tions should also include this effect. We take into ac- 
count the Coulomb interaction by using the exact form 
of the two-particle wave-function 12ofl . In the procedure 
of including the Coulomb effects we should accordingly 
modify the form of the reference function I|27J1 that is 
used to extract the HBT radii. Since there is no analytic 



9 




q [GeV/c] 

long L 




0.08 0.1 
q [GeV/c] 

long L 



FIG. 8: An example of the projected pion correlation func- 
tions for the Blast-Wave model with resonances, a — —0.5. 
Model parameters are the same as for Fig. The correlation 
functions include pion pairs with transverse-momentum in the 
range: 0.25 GeV < kr < 0.35 GeV. We show the projections 
of the correlation function (symbols) and the projections of 
the 3-dimensional fit (lines). The top plot shows projections 
on the (/out axis, i.e., the 3-dimensional correlation function 
has been integrated in two other directions over some range. 
The middle plot shows the projection on g s id c , and the bot- 
tom plot the projection on qi ons - Circles (solid lines) show the 
function (fit) integrated in two other directions in the range: 
\qi\ < 2.5 MeV, triangles (dashed lines) for \qi\ < 12.5 MeV, 
and squares (dotted line) for < 32.5 MeV, where i = out, 
side or long. 



formula which parametrizes the Coulomb effects exactly, 
we follow the procedure of CERES [12], STAR HH and 
PHENIX ^3], which is based on the following assump- 
tions: the Coulomb interaction and the wave-function 



FIG. 9: The same as Fig. [&] for the analysis including the 
Coulomb interaction. 

symmetrization factorize and, moreover, the Coulomb in- 
teraction part of the function can be replaced by the av- 
eraged Coulomb wave-function. With these assumptions 
the Coulomb interaction may be separated from the in- 
tegration in the numerator of lt25]l and one obtains the 
Bowler-Sinyukov formula 0, l4fil| 

C(q, fe) = (1 — A) + AX cou i(g inv ) [l + cxp (-R 2 out ql ut 

~ ^sidc9sidc ~ ^long<Zlong)] ! (28) 

where K cou i(qi nv ) is the squared Coulomb wave function 
integrated over a static gaussian source. We use, follow- 
ing the STAR procedure ^l|]> the static gaussian source 
characterized by the widths of 5 fm in all three direc- 
tions. The 3-dimensional correlation function with the 
exact treatment of the Coulomb interaction, calculated 
according to Eqs. 11411 and 125L is then fitted with this 
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FIG. 10: Ratio -R ou t/-R s idc for the Cracow model with T = 165 
MeV (circles) and the Blast- Wave model with resonances and 
a = —0.5 (squares). The STAR data are represented by open 
circles. 

approximate formula and the HBT radii are obtained. 
They can be compared directly to the experimental radii. 
We note that this way of comparing experimental data 
and theoretical predictions is very reasonable in the sense 
that the same experimental and theoretical observables 
are compared. The HBT radii may be obtained also from 
the fit to the correlation function 1270 , which provides the 
way for the experiments to judge the systematic uncer- 
tainty of determining the results by using the Bowler- 
Sinyukov procedure. 

IV. RESULTS 

The parameters of each model are fixed by fitting the 
single-particle p_L-spectra of pions and kaons to the ex- 
perimental data. An example of such a fit is shown in 
Fig. OH The values of the parameters are listed in the 
captions of Figs. 14171 Thus our analysis of the HBT cor- 
relations has no extra parametric freedom. 

Our basic results for the pion interferometry are 
shown in Figs. 141 71 They were obtained by the pro- 
cedure described in detail in Sect. IIII1 i.e., by fitting 
the 3-dimensional two-particle correlation functions. In 
Figs. 14171 we show the intercept A and the HBT radii 
-Rout, -Rsidc, and i?iong as functions of the transverse mo- 
mentum of the pion pair. The squares correspond to 
the full calculation with resonances, the down-triangles 
show the results obtained in the calculation without res- 
onances, the up-triangles show the results obtained with 
resonances and with the Coulomb-aware fit made accord- 
ing to the Bowler-Sinyukov formalism while the 
circles show the data of the STAR collaboration from 
Ref. j^J. The first immediate observation is that the in- 
clusion of resonances increases the radii by about 1 fm. 
This is expected, since the resonances travel some dis- 
tance from their place of birth on the freeze-out hyper- 



surface before they decay into pions. The typical scale is 
set by the resonance life-time which is about 1 fm/c. 

Browsing through Figs. 1 1171 we note that the agree- 
ment with the data changes with the selected model of 
expansion. Out of the four models tested, by far the 
best results are obtained for the Blast- Wave model with 
resonances and a = —0.5. This shows that the hypersur- 
faces with p decreasing with time (c/. Fig.^l are favored. 
Note that the iiiong radius is particularly sensitive to the 
a parameter, with a — —0.5 giving the right result, while 
increasing a spoils the agreement. The model values of 
the intercept A shown in Figs. 14171 are too large compared 
to the data, which simply reflects the fact that we do not 
take into account the effect of secondary pions coming 
from the weak decays, as well as the contamination of 
the pion sample by misidentified particles in the experi- 
ment. 

The Cracow model and the Blast- Wave model with 
a = 0.5 have very close predictions, as expected from the 
similarity of the hypersurfaces, cf. Fig. We note that 
in all considered models the qualitative behavior of the 
dependence of the radii on kx is correct. The Coulomb 
corrections evaluated with the Bowler-Sinyukov formal- 
ism are small, of the order of a small fraction of a fermi. 

Our method of determining the HBT radii from 
the model involves the calculation of the complete 3- 
dimensional correlation function. First we use the free 
wave function of Eq. 12 -ill in Eq. 12211 and compare it to 
the gaussian fit made according to Eq. 127fl . The results 
are shown in Fig. where we plot the projections of the 
correlation function itself, as well as the projections of the 
3-dimensional fit. The deviations between the function 
(symbols) and the fit (curves) reflect the fact that the un- 
derlying two-particle distributions are not gaussian and 
produce a non-gaussian correlation function. One can 
also see that the increase of the integration regions in 
the complementary directions leads to better agreement. 
Since this is an important finding, we restate this obser- 
vation: In a fixed kr bin we take the 3-dimensional func- 
tion. We choose one of the directions, say g ou t, and inte- 
grate over the remaining two directions, q S ido and Qlong, 
within the specified range. When the range is very nar- 
row, this corresponds to slicing the 3-dimensional func- 
tion along the line q sl d c = 9iong = 0. This gives the circles 
in Fig. Next we repeat the same projection prescrip- 
tion but for the gaussian fit to the correlation function, 
which results in the solid lines. The lines deviate from 
the circles within a few percent showing that the gaus- 
sian approximation works at that level. Increasing the 
integration range in the complementary directions results 
in a much better agreement, which is manifested by the 
overlap of the squares to the dotted lines. One can also 
see that even though the detailed shape of the correlation 
function is not exactly reproduced when the integration 
region is smaller (circles and triangles), the overall width, 
and hence the radius, is described well. 

Now we come to the analysis of the correlation function 
taking into account the Coulomb interactions. Explicitly, 
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FIG. 11: The separation distributions for pion pairs from the Blast- Wave model with resonances and a = —0.5 (black circles). 
On the left-hand-side the plots for pairs of primordial pions are shown. On the right-hand-side, the plots for all pions are 
shown. The lines show the separation distribution which is the result of the fitting of the corresponding correlation function 
by a gaussian parameterization. 



we use the Coulomb wave function of Eq. (12ofl in Eq. J22I) 
and compare it to the Bowler-Sinyukov formula l(28l) . We 
perform the same procedure as above and the results are 
shown in Fig. El We observe that the Coulomb interac- 
tions dig holes at low values of q, which is the well-known 
result of the long-range repulsion. 

In Fig. El we show the ratio i? ou t/-Rside for several 
models. The very good agreement with the STAR data 
(open circles) is obtained for the Blast- Wave model with 
resonances and a = —0.5. This behavior is already ex- 
pected from the inspection of Figs. 14171 The Cracow 
model with T = 165 MeV describes well the kr depen- 



dence of the experimental ratio, giving the magnitude of 
the ratio smaller by about 15%. 



V. EFFECTS OF RESONANCES IN THE 
CORRELATION FUNCTION 

Figure El shows the separation distributions for the 
Blast- Wave model with a = —0.5 and for the bin 0.25 
GeV < k T < 0.35 GeV. On the left-hand-side the dis- 
tributions of pairs constructed only from the primordial 
pions are shown. It can be seen that these distributions 
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FIG. 12: Space distribution of the produced pion pairs. The 
solid lines show separation distribution for all pairs, circles 
show the pairs of primordial pions, up-triangles show the pairs 
which contain one pion from the p decay, but no pion from 
the u> decay, squares show the pairs where one pion comes 
from the U) decay, finally down-triangles show the pairs where 
one pion comes from decays of resonances other than p and 10 
and the other from this group plus the primordial pions, see 
text for a more detailed description. Blast- Wave model with 
resonances, a = -0.5, 0.25 GeV < k T < 0.35 GeV. 



are cut off at the value determined by the p m ax parameter 
of the Blast- Wave model, since by definition no pair can 
have a separation greater than 2p max in the out and side 
directions. The correlation function is the Fourier trans- 
form of these distributions. 



As described in Sect. IIII PI 
the correlation function may be fitted with the gaussian 



formula (23 and the result of the fit can be used to find 
the presumed gaussian distribution l(2fijl . The separation 
distributions obtained with the help of such a procedure, 
i.e., from the gaussian fits to the correlation function, are 
shown as thin lines in Fig. El We note that the gaus- 
sian approximation works reasonably well at low values 
of the separation radii with some deviations in the case 
of long direction at large values of n ong . In that case a 
tail results from the chosen model of expansion. In our 
boost-invariant model the distribution extends to infin- 
ity, however, the observed drop is caused by the presence 
of the homogeneity length in the system. 

On the right-hand-side of Fig. El we show the cor- 
responding plots including all pairs of the pions created 
before 500 fm/c, i.e., the primordial pions and other pi- 
ons coming from almost all strongly-decaying resonances. 
The effect of the resonances is very clearly visible in a 
long-range tail in the out and side directions. One can 
see how the model curve departs from the gaussian fit 
around 17 fm. For the long direction the tail is also pro- 
duced and the effect adds up to the tail produced by the 
expansion model. Thus, it is clear that the gaussian fit 
to the correlation functions has no way to describe the 
long-range tail of the distributions. The tails have, how- 
ever, an effect on the correlation function and show up 
as peak at low q mv which is seen in FiglHl as the differ- 
ence between the data (circles) and fit (solid lines). This 
results in the lowering of the A parameter of the gaussian 
fit Q. 

It is interesting to study the long-range tails of the sep- 
aration distributions in more detail. Figure 1X21 shows the 
anatomy of the separation distributions divided into sev- 
eral components. The pions are divided into four groups: 
primary, those coming from p decays, those coming from 
u> decays, and other coming from decays of other reso- 
nances than p or lu. In the plot we show the distributions 
of pairs constructed from pions belonging to the classes 
defined above. First we consider pairs where both pions 
are primary (circles). One can see that in the case of out 
and side directions these pions are concentrated near the 
origin with 2p max ss 18 fm providing the cutoff. There 
is no such cutoff in the long direction where we can see 
the falloff resulting from the homogeneity length. Next 
we consider the pairs where one of the pions comes from 
the p decays and the other from primary pions or pi- 
ons coming from decays of resonances other than u> and 
p (up-triangles, labeled as "primary, other - p"). These 
pairs are responsible for the increase of the strength of the 
source in all three directions. In the out and side direc- 
tions, they cause the swelling of the source. The curves 
corresponding to the pairs where one of the pion belongs 
to the group "other" and the second to the groups "other, 
primary" (down-triangles) show a very similar behavior 
to the "primary, other, p" case. Finally, we show the 
pairs where one of the pions comes from the u> decay 
(squares). In all three directions we observe the long- 
range tails caused by the long lifetime. These pairs pro- 
duce a non-gaussian character of the correlation function. 
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PIG. 13: The 1-dimensional correlation function plotted as a function of qi m for the Blast-Wave model with resonances, 
a = —0.5, 0.25 GeV < kr < 0.35 GeV (open circles). Panel a) gives the "primary-primary" case, b) - "primary, other - p", 
c) - "primary, p, other - lo" , and d) - "other, primary - other". For comparison dashed lines show the best fit of the gaussian 
formula, with the corresponding values of i?i n v 



It can also be seen, that the total distribution (solid lines) 
can be well approximated by a combination of a gaussian- 
like core at low r and the long-range non-gaussian halo. 

In order to illustrate the influence of the resonances on 
the correlation function itself, in Fig. El we show how 
different types of pairs contribute to the observed corre- 
lation as a function of q- mv . It is clearly seen that none of 
the contributions is well-reproduced by a gaussian. This 
feature is most prominently seen for the pairs which con- 
tain at least one pion from the oj decay: in this case 
the correlation function has a peak at low qi nv , which 
is expected as the spacetime distribution for these pairs, 
seen in Fig. El is exponential at large r. Even though 
the gaussian fit fails to describe the detailed behaviour of 
the functions, especially at low values of q mv , it serves as 
a tool for extraction of the HBT radii. One can see, that 
the smallest radius is obtained for the primordial pairs, 
as expected. However, they provide only about 10% of 
the correlation effect. Pairs, which contain a pion from 
any strongly-decaying resonance, except for p or u, show 
a size larger by approximately 1 fm. They account for at 
least 30% of the correlation effect. Pairs containing the p- 
decay product show slightly larger increase of the radius 



(also about 1 fm) and provide the largest contribution to 
the correlation effect, about 40%. Pairs containing the 
w-decay products give the largest size, as expected, but 
their contribution to the correlation function is sharply 
peaked at small qinv, which results in the decrease of the 
A parameter, but does not influence the width of the cor- 
relation function (and therefore the obtained radius). 



VI. PREDICTIONS FOR THE KAON 

In this section we show our predictions for the kaon 
HBT radii from the Blast- Wave model with resonances. 
The results obtained with the free wave function l|23jl 
and formula 12211 are shown in Fig. El We also give the 
results for the pions and notice that the kaons nicely ex- 
tend the curves to higher values of the transverse mass. 
The results are plotted as a function of tot of the pair. 
This choice follows from the expectation that the radii 
are influenced mostly by the flow, which results in the 
approximate tot scaling 03| ■ The error bars indicate 
the errors of the gaussian fit to the THERMINAT0R simula- 
tion. We note that for all radii the scaling holds within 
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FIG. 14: Predictions for the kaon femtoscopy in the Blast- 
Wave model with a = —0.5. The model parameters are 
T = 165.6 MeV, fi B = 28.5 MeV, r = 8.55 fm, p max = 8.92 fm, 
vt = 0.311 c. The first four points from the left are predic- 
tions for pions, the three points on the right - for kaons. 



the chosen model of freeze-out as well as the detailed 
influence of resonance decays. Concerning the choice of 
the expansion model we have found that the freeze-out 
geometry where the transverse size decreases with time 
works by far the best, allowing for a uniform description 
of the pr spectra and the HBT correlation radii. Such 
a model has the features similar to those found in the 
advanced hydro calculations, however in our study the 
resulting time and size scales are shorter. We have found 
that the i£i on g radius is particularly sensitive to the de- 
tails of expansion, which helps to discriminate between 
various cases. Our calculations were done for the Cra- 
cow single-freeze-out model as well as for the Blast- Wave 
model including all resonance decays, and with a modi- 
fied shape of the p — t freeze-out curve. We have achieved 
a satisfactory and uniform agreement for the description 
of the data, as can be judged from Figs. 3 and 7. 

Our calculations used the code THERMINAT0R, which 
is a Monte-Carlo implementation of the hydro-inspired 
models with single freeze-out. The use of the Monte- 
Carlo technique allowed us to study in a greater detail the 
effects of resonances on the shape of the correlation func- 
tions. Since our freeze-out temperature is large, we need 
to include practically all resonances, as in the studies of 
particle abundances and momentum spectra. We have 
explicitly found non-gaussian features of the pion cor- 
relation functions caused by the long-living resonances, 
mainly the u> meson, confirming earlier expectations. In 
addition, we have carefully discussed their quantitative 
role for the extraction and interpretation of the HBT 
radii, as well as the shape of the correlation functions and 
the separation distributions. In short, we hope that our 
analysis provides a very useful "vivisection" of the pion 
HBT problem, helpful in the understanding of the un- 
derlying space-time picture of relativistic heavy-ion col- 
lisions. We find that the pion HBT data from RHIC are 
fully compatible with the single freeze-out scenario. 

Finally, we give predictions for the HBT radii of kaons, 
which should be measured shortly. The results for the 
kaons exhibit the ray-scaling proposed in Ref. (IfiL l47l| . 



APPENDIX A: CORRELATION FUNCTIONS 
FOR PRIMORDIAL PARTICLES 



these uncertainties. The calculation presented in Fig. 1141 
uses the model which was most successful for the pions. 
Note that all parameters were fixed by the pionic sec- 
tor and the predictions for kaons involve no parametric 
freedom. 



In the case of the primordial particles one may obtain 
analytic expressions for the correlation functions which 
are very useful as the reference point for the Monte-Carlo 
method. They can also be used for testing numerical 
calculations. 



VII. CONCLUSIONS 

In this paper, a class of hydro-inspired models with 
single freeze-out was used to analyze in detail the pion 
correlation functions, with an emphasis on the role of 



1. Side correlation function 

In the case of the side correlation function, one may 
choose the direction of the average transverse momentum 
to be parallel to the r x - axis and the difference of the two 
momenta to be parallel to the r y - axis. In this reference 
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frame, at zero rapidity (y = 0) we have 

k = (kj_,0,0), q = (0, q sidc , 0), 

and the phase of the Fourier transform of the emission 
function has the form 

i q ■ x = -i q sldc y = -i g S ido P sin (j>. 

Performing similar steps as those which led us to Eq. <@J 
we obtain the Fourier transform of the emission function 
in the following form {q = q s i de ) 

oo 1 

S side (k ± ,q) = ^EW n+1 e"^| d( p(C)f(C) 



n = l 



dp 



x <^m±— K\ [n/3r7i^cosha^]Xo [n(3k±smha±, qp{()] 
df 

—p± — K a [n/3m±cosha±] X\ [n/3/c^sinhax, qp(()] 
at, 



For 6 = they are reduced to the Bessel functions Io(a) 
and h(a). 



2. Out correlation function 



In this case we may choose the coordinate system 
where 



fc = (fcL,0,0), q= (^ut, 0,0), 



and the phase of the Fourier transform is 



i q ■ x = iq°t~i q sidc x 



i q f coshan — i q sldc p cos < 



Here we introduced two functions: 

2tt 



If 

Xo(a,b) — — / dip exp (a cos tfi — ib sin ( 

2-7T J 



The energy difference q° in Eq. (I A lfi is obtained from 
the formula 



q° = \Jm 2 + (fc_L + <?sidc /2) 2 - ^m 2 + (fc ± - <7sidc /2) 2 



and 



Z7T 

Xi(a,b) = — / dtp cost/) exp(acos0 — iftsinc 

27T 7 



and the Fourier transform of the emission function reads 

(q = <7side) 



S ou t(k±,q) 



oo 1 
n— i n 



-p±-rp K [n[3m±cosha± — f] 7i {nf3k±sinha± — iqp] 
dQ 



3. Long correlation function 



and (q = giong) 



Similarly to the two previous cases, we find 



fe = (fe_L,0,0), q= (0,0,Qi on g), 



i 5 • x = — i giong z = — i gi ong fsinhan , 



00 /• 

W*x,?) = ^E(t)" +1 ^ / d( p(C)f(C) 

n=l 

x I m±-JpfCi [n/3mj_cosha±, gf] 1q [n/3fcj_sinha!_L] 

df I 
— p±— ICo [n{3m±cosha±, qf } 1\ [n/3fc^sinhajj > . 
d( J 
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The functions /Co and 1C\ are defined by the integrals: For 6 = these functions reduce to the modified Bessel 

functions K and K\. 



/Co(a,6) = i J dacxp(— acosha — ibsmha) 



and 



/Ci(a,6) = - j da cosha exp(— acosha — ibsinha). 
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